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A FORTRAN  program  for  linear  predictive  spectral  analysis  of  univariate 
complex  data  with  disjoint  equi-sized  pieces  is  presented.  This  result  has 
application  to  data  taken  with  a device  that  periodically  goes  off-line,  and 
to  maximum  entropy  beamformlng  on  a limited  number  of  equi-spaced  linear 
array  elements.^ 
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FORTRAN  PROGRAM  FOR  LINEAR  PREDICTIVE  SPECTRAL  ANALYSIS  OF 
UNIVARIATE  COMPLEX  DATA  WITH  DISJOINT  EQUI-SIZED  PIECES 

INTRODUCTION 


In  reference  1,  the  theoretical  basis  for  estimating  spectra  via 
linear  predictive  techniques  from  univariate  complex  data  with  bad  data 
points  was  presented.  Also,  a program  was  given  there  for  spectral 
analysis  of  real  data  with  gaps.  Later,  in  reference  2,  a program  for 
complex  data  with  no  bad  data  points  was  presented.  However,  neither 
program  is  capable  of  handling  complex  data  with  gaps.  This  present 
document  remedies  this  situation  for  the  case  where  the  disjoint  pieces 
of  complex  data  each  have  the  same  number  of  points. 

Complex  data  can  arise  when  a real  process  is  complex-demodulated 
to  zero  frequency  and  sampled  at  a low  rate  for  purposes  of  ease  of 
processing  and  computation.  Alternatively,  if  a real  process  is  sub- 
jected to  Fourier  transformation  into  the  frequency  domain,  complex 
coefficients  result.  Both  of  these  examples  are  of  frequent  occurrences 
in  practical  applications. 

Data  segments  of  equal  length  can  occur,  for  example,  when  a 
recording  device  is  periodically  taken  off-line  (for  calibration  purposes, 
perhaps),  or  if  large  periodic  bursts  of  noise  occur.  An  important 
frequency-domain  application  occurs  for  an  equi-spaced  line  array  of 
limited  extent.  Then  the  number  of  elements  in  the  array  is  the  maximum 
size  of  each  disjoint  piece  of  frequency  coefficients.  If,  for  a partic- 
ular time  segment  of  array  outputs,  we  look  at  the  Fourier  coefficients 
at  one  frequency,  we  can  attempt  linear  prediction  in  space;  averaging 
the  prediction  error  over  time  then  yields  the  extrapolation  effect 
inherent  in  maximum  entropy  processing  (see  reference  3).  This  will 
result  in  the  array  appearing  to  have  a longer  length  than  its  true 
length,  and  thereby  yield  improved  angular  resolution  of  plane-wave 
arrivals.  This  approach  is  analogous  to  the  improved  frequency  resolu- 
tion obtainable  from  limited  time  data. 

The  input  parameters  to  the  program  are  fully  explained  in  the 
comment  statements.  For  a single  piece  of  data,  the  number  of  disjoint 
pieces,  ND,  can  be  set  equal  to  1.  A sample  run  is  presented  after  the 
program.  For  application  to  linear  predictive  beamforming,  delete  the 
statement 


X(I,  J)  = X(I,  J)  - AVE 


in  loop  3 of  SUBROUTINE  CXDISJ. 
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LINEAh  predictive  SPECTRaL  ANALYSIS  FOR  UNIVARIATE 
COMPLEX  DATA  WITH  DISJOINT  EQUI-SIZED  PIECES 
C EQUATION  REFERENCES  ARE  TO  NUTTALL»  NUSC  TR  5303»  26  MARCH  1976 
C user:  CHANGE  LINES  26  ANq  40»  aND  REPLACE  SUBROUTINE  DATA 
C rj  = NUMBER  OF  COMPLEX  CAtA  POImTS  IN  EACH  PIECE  I INTEGER  INPUT 
C lJU  = NUMBER  OF  DISJOINT  pIECESI  INTEGER  INPUT 
C X(1,1),,.X(N»1),...»X(1»nD)...X(N»MD)  = COMPLEX  INPUT  DATAI 
C ALTERED  ON  OUTPUT 

C PMAX  = maximum  ORDER  OF  pIlTERi  PMAX.LT.NJ  INTEGER  INPUT 
C UF  = SI/E  OF  FFT  (MUST  0E  A POWER  OF  2 TO  USE  MKLFFT)  > INTEGER  Ik, PUT 
C AVE  = complex  sample  MEArt  oF  INPUT  DaTAI  OUTPUT 
L Pn  = SAMPLE  POWER  OF  lUPtjT  DATAI  OUTPUT 

C AIC  = AKAIKE’S  INFORMATION  CRITERIONi  OUTPUT 
C PBEST  = BEST  ORDER  OF  FIlTER*  INTEGER  OUTPUT 
C SPBFsr  = (EO  H-7)/MF  FOR  prPBESTI  OUTPUT 
C SPMAX  = (EQ  H-7)/NF  FOR  PsPMAXI  OUTPUT 

C A{1) »..,»A(PBEST)  = complex  PREDICTIVE  FILTER  COEFFICIENTS  = 

C AdIPBEST A (PEEsTiPBEST)/  OUTPUT 

c RHO(i) »,..»rho(pmax)  = Complex  normalized  correlations:  output 
c xx(i) »,..»xx(MF)  = Scaled  spectrum,  whose  sum  = sample  power:  output 

C C0Sl(l)».,.»C0SI(NF/4+l)  = QUARTER  COSINE  TABLE  FOR  FFT  PURPOSES 
C Y(N,Nb)  IS  A REQUIRED  COMPLEX  AUXILIARY  ARRAY 
C YY(r'F)  IS  A REQUIRED  AUXjLlARY  ARRAY 

C ON  OUTPUT,  X(l,l) ,,,,,X(pMAX,l)  = A ( 1 I PMAX) » . . . » A (PMAX IPmAX ) 

C ON  OUTPUT,  Y(l,l) ,.,,,Y(pMAX,l)  = A ( 1 1 1 ) » . . . » A ( PMAX I PMAX ) 

PARAMETER  N = 1C,  r;D  = 20,  PMAX  r 6,  NF  r 512,  MF41=MF/4  + 1 

INTEGER  PBEST»P 
real  SPBEST,T1,T2 

complex  X(N,ND) ,Y(N,Nr) jACPmAX) , RHO (PMAX) , AVE, G»T 
DIMENSION  XX(NF) ,YY(Nf) ,C0SI (NF41) , AIC (PMAX) ,AIC0 (2) 

EQUIVALENCE  ( AIC ( 1 ) , AjCO (2) ) 

C PRINT  OUT  VALUES  OF  PARAMETERS 
I=N 
U=ND 
K=PMAX 
L=NF 
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PRINT  If  IfJfKfL 

I FORMAT(lHlf»  N =• f I6» lOX» *ND  =♦ f I6f lOXf ’PMAX  S* f I4f lOXr »NF  =»,IS) 
C COMPLEX  INPUT  data  IN  X ( 1 f 1 ) • • .X (Nf 1 ) f , , , » X ( 1 f NO) . • .X (Nf ND) 

CALL  data 
PRINT  2 

k FOKMAT(/»  COMPLEX  INPuT  DATA!*) 

DO  3 J=lfND 
PRINT  34f  J 

34  F0KMAT(»  piece  NUMBFR»fl3) 

PRINT  4f  (XdfJlf  1=1, N) 

F0RMAT(4(Eie.8,E15.P) ) 

EVALUATE  FILTER  COEFFICIENTS 
■ ’call  CXDISJ 
PRINT  5f  AVE 

t)  FORMAT(/»  COMPLEX  MEAN  OF  INPUT  DATA  = ( » »E13.8f » f ♦ f El3.Bf  • ) • ) 
PRINT  6f  PO 

o FORMAT(/»  sample  POVEr  oF  INPUT  DATA  =»fE13,8) 

PRINT  7 

7 FORMATC/'  AKAlKE  INFORMATION  CRITERION: ♦ /9X f »P* f IIX ,♦ AlC (P) • ) 

PRINT  8f  (P,AlC(P)f  P=0,IA) 
t>  FORMAT(I10,E20,8) 

PRINT  9f  PBEST 
9 FORMAT(/*  PBEST  =*fl3) 

PRINT  10 

1''  FORMAT(/»  PARTIAL  CORRELATION  COEFFICIENTS:*/ 

$9Xf *P»f9Xf *RE  A(PIP) »,7X**IM  A(Plp)») 

PRINT  Ilf  (PfY(Pfl)f  P=lflA) 

II  FOKMAT(I10fE20.8fE16.8) 

IF(P8EST,EQ,0)  GO  TO  i2 
PRINT  13 

13  format:/*  predictive  filter  coefficients  for  pbest:*/ 

sgXf  ♦KSTXf  *RE  A(KlPPEsT)*f3X»*lM  A(KlPBEST)*) 


PRINT  Ilf  (PfA(P),  P= 

If PBEST) 

niN 

,1.  ■ 

. 1 

3 


12  PRINT  14 

14  format(/»  normalized  correlations:*/ 

$7X» 'DELAY*, 9X» ‘RE  RHO*»lOX,*IM  RHO*) 

P=0 

T=(1.»0.) 

PRINT  ll»  P»T 

PRINT  ll»  (P,RHO(P),  P=1»IA) 

C EVALUATE  SPECTRAL  ESTIMAjE 
CALL  SPECT 
PRINT  15 

15  format:/'  power  spectrum»  starting  from  zero  frequency  {Bind:*) 

PRINT  16»  (XX(I)»  1=1, NF) 

16  F0RMAT(2X»10E13.S) 

PRINT  17 » SUM 

17  format:/*  sum  of  spectrum  values  =*,E13,8) 

PRINT  18»  PO 

1ft  format:*  sample  power  of  input  =*,E13,8) 

c 

SUBROUTINE  DATA 

C THIS  subroutine  GENERATES  COMPuEX  DATA 

define  IRAND=I*5**15-I-  ( { I-SIGN:  1 1 1*5**15) ) /2)  ♦34359738367 
define  RANDsFLOAT: I )/34359738367. 

1=5281 
6=:. 65,. 65) 

DO  1 J=1»ND 
T=:0.rO.) 

DO  2 K=l»200  B will  DISCARD  THESE  INITIAL  POINTS 

I=IRAND 

Tl=RAND-.5 

I=1RAND 

T2SRAN0-.5 

2 T=G^T+CMPLX:T1»T2) 

DO  3 K=1»N 
I=1RAN0 
T1SRAND-.5 
I=IRAND 
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PY 


E 


T2=RAND-.5 
T=e*T+CMPLX(Tl»T2) 
i X{K»J)=T 

1 CONTINUE 

RETURN 
C 

SUBROUTINE  CXDISJ 

C THIS  SUBROUTINE  COMPUTES  AlC»  PBEST»  ALL  THE  FILTER 
C C0EFFICIENT5»  AND  THE  FORMALIZED  CORRELATIONS 
rOOBLE  PRECISION  SAR»SAI»SB 
I=N 

K=PMAX 

IA=N-1 

IF(PMAX,GE.N)  PRINT  1,  K»I»1A 

1 FORMAT(/'  PMAX  =»,!*<»»  IS  TOO  LARsE  FOR  NUMBER  OF  DATA** 

•i*  points  N =«»I5»*»  search  limited  TO  P =*»I4) 

IA=MIN(IA»PMAX) 

rACr4./(N*ND)  Ci  FAc=0.  WOULD  FORCE  PBEST  EQUAL  TO  I A 
C COMPUTE  SAMPLE  MEAN 
A\/E=(C«  »0*  ) 

DO  2 J=l*ND 
DO  2 I=1*N 

2 A»/L=AVE+X(I,J) 

A W E=CMPLX ( re AL ( a VE ) / ( N*ND ) » a IMAG ( a VE ) / ( N*ND ) ) 

C SUBTRACT  SAMPLE  MEAN  AND  COMPUTE  SAMPLE  POWER 
P0=0. 

DO  3 J=l»ND 
DO  3 1=1 »N 

X(I»J)=X(I»J)-AVE  C to  keep  sample  mean,  delete  this  CARD 
Y(I,J)=X(I,J) 

3 P0=P0-*-REAL  ( X ( I , J ) ) **24>AIMA6  ( X ( I * J ) ) *«2 
P0=P0/(N*ND) 

C BEGIN  RECURSION 

AIC(0)=LOG(PO) 

AICMIN=AIC(0) 

P8EST=0 
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SPMAX=PO/NF 
SP0EST=SPMAX 
DO  4 P=1»IA 

C CALCULATE  CROSS-GAINI  EQ  155 
SARrO.on 
SAIsO.DO 

sa=o,Do 

L=P+1 

DO  5 J=1»ND 
DO  5 I=L»N 
n=REAL(X(I,J) ) 

T2=AIMAG(X(I,J) ) 

T3=REAL(Y(I-1»J) ) 

T4=AIMAG(Y(I-1,J) ) 

SAR=SAR+T1*T3+T2*T4 
SAI=SAI4-T2*T3-T1*T4 
b SB=SB+T1**2+T2*»2+T3**24T4**2 

T1=2,*SAR/5B 
T2=2.*SAI/SB 
G=CMPLX{71»T2) 

C CALCULATE  FILTER  COEFFICiEnTSi  EQS  160«148.  STORE  IN  X ( 1 » 1 ) . , ,X (p » 1 ) 
X(P*1)=G 

IFlP.EQ.l)  GO,  TO  6 
L=P/2 

DO  7 I=1»L 

T=X(I»1)-G*C0NJG(X(P-I»1) ) 

X ( P- I » 1 ) =X ( P-I , 1 ) -G*CoN JG ( X ( I » 1 ) ) 

X(I»1)=T 

CALCULATE  NORMALIZED  CORRELATION!  EQ  149 
T=X(P»1) 

IF(P,EQ,1)  GO  TO  8 
L=P-1 

CO  9 1-ltL 
T=T-»-X(I»l)*RHO(P-n 
PH0(P)=T 
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c calculate  akaike's  information  criterion 

T1=1.D0-4.D0*(SAR**2+SAI**2)/SB**2 
SPMAX=SPMAX*T1 
AIC(P)=L0G(SPMAX)+FAC*P 
IF(AIC(P) .GE.AICMIN)  GO  TO  lO 
AICMIN=AIC(P) 

PbEST=P 
SPbEST=SPMAX 
DO  11  I=1»P 

11  A(I)=X(I»1) 

10  IF(P,EQ,IA)  60  TO  12 

C UPDATE  FORWARD  AND  BACKWARD  SEQUENCES I EQ  153 

L=P+1 

DO  13  J=1»ND 
DO  13  I=N»L»-1 
T=X(I»J)-G*Y(I-1»J) 

Y(I»J)=Y(I-1, J)-CONJG(G)*X{I,J) 

13  X(I*J)=T 

4 YCP»1)=G 

12  Y(IA»1)=G 
IF(PBEST.EQ,IA)  GO  TO  14 

C COMPUTE  extrapolated  NORMALIZED  CORRELATION 
C coefficients  from  PBEST+1  to  PMAXI  Eq  165 
L=PBEST+1 
DO  15  P=L»IA 
A(P)=(0.»0.) 

T=(0,»0.) 

DO  16  I=1»PBEST 
16  TsT-fACD^RHOCP-I) 

15  RHO(P)=T 

14  RETURN 
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SUBROUTINE  SPECT 

C THIS  subroutine  COMPUTES  THE  POWER  SPECTRUM  FOR  PBESTI  IT  IS  SCALED 

c SUCH  that  the  sum  of  values  computed  should  equal  the  sample  power 

XX(1)=1. 

•Yr(i)=o. 

IF(PBEST.EQ,0)  GO  TO  i 
DO  2 I=1»PBEST 
XX(I+1)=-REAL(A(I) ) 

2 YY(I+1)=-AIMAG(A(I) ) 

1 L=P0EST+2 

DO  3 I=L»NF 
XX{I)=0. 

3 YY(I)=0. 

CALL  QTRCOS(COSI»NF) 

L=1.4427*L0G(NF)+,5  Q LOg2(NF) 

CALL  MKLFFT(XX,YY,C0Si»L»“1) 

SUM=0. 

DO  4 I=1»NF 

XX  ( I ) =SPBEST/  UX  ( I > ♦♦a+YY  1 1 ) **2 ) 

4 SUM=SUM+XX(I) 

RETURN 

END 


SUBROUTINE  QTRCOS(C»N) 
DIMENSION  C(l) 
N41SN/4+1 
SCLS6.283165307/N 
DO  1 I=1»N41 

1 C(I)=COS((I-l)*SCL) 
RETURN 
END 
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SUtiROUTINE  MKLFFT(X»Y,CC»M,ISN) 

Ci,^iENSION  X(l)  »Y{1)  »Cc(l)  fLda) 

F^UlVALEf'JCE  (L12»L(  J ) ) ► {L11,L  (2)  ) » (L1P»L{3)  ) » {L9,L(4)  ) r (LBr L(5)  ) » 
1<l7»L(6) ) » (L6»L(7) ) » (I  5,L(8) ) » {L4,L(9)) » (L3»L(10) ) , {L2fL(ll) ) , 
2(l1»L(12) ) 

^ L)H=rxj/4 
f J4P1=ND4+1 
I.U4P2=ND4P1  + 1. 
r J2P2=ND4+hD4P2 
DO  B LO=l»M 
L.»1X=2**(M-L0) 

LlX=2*L^’X 

ISCL=N/LIX 

LO  B LM=1»LMX 

IAkG=(LM-1)*ISCL-*-1 

IF(IARG.LE.ND4P1)  GO  yO  4 

C=-CC{ND2Pa-IARG) 

S=ISN*CC(lARG-ND4) 

GO  TO  6 
*+  C=CC(IARG) 

S=i5N*CC (ND4P2-IARG) 

6 DO  8 LI=LIX»^J»LIX 
J1=LI-LIX+LM 
J2=J1+LVX 
T1=X( J1)-X(J2) 

T2=Y(J1)-Y{J2) 

X( J1)=X(J1)+X(J2) 

Y(J1)=Y(J1)+Y(J2) 

X{ J2)=C*T1-S*T2 
Y(J2)=C*T2+S*T1 

6 continue 
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DO  40  J=l»12 
L(J)=1 

IF(J-M)  31,31»40 
31  L(J)=2**(M+1-J) 

40  CONTINUE 
JN=1 

DO  60  J1=1»L1 
DO  60  J2=Jl»L2rLl 
DO  60  J3=J2fL3»L2 
DO  60  J4=J3tL‘*rL3 
DO  60  J5=J4»L5»L4 
DO  60  J6=J5»L6»L5 
DO  60  J7=J6»L7»L6 
DO  60  J8=J7»LB,L7 
DO  60  J9=J8,l9»L8 
DO  60  J10=J9»L10»L9 
DO  60  J11=J10»L11»L10 
DO  60  JR=J11»L12»L11 
IF(JN-JR)  51, 51 » 52 

51  R=X{JN) 

X(JN)=X( JR) 

X(JR)=R 

FI=Y(JN) 

Y{JN)=Y( JR) 

Y(JR)=FI 

52  JN=JN+1 
60  CONTINUE 

RETURN 

END 
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COMPLEX  MEAN  OF  INPUT  DATA 


(-.66319221-01 40282205-01) 


bAMPLt  POWER  OF  INPUT  DATA  = . 1224ai06<»-0l 

AKrtIKE  information  CRITERION: 

P AIC(P) 

0 ,202786204-00 

1 -,601698104>01 

2 -,800728504-01 

3 -,798759594-01 

4 -,797029574-01 

b -,796545034-01 

b -,795190504-01 

PBEST  = 1 

PARTIAL  correlation  COEFFICIENTS: 

P RE  A(P»P)  IM  A(P»P) 

1 ,656078984-00  .659105584-00 

2 ,71265785-01  ,719l8564-0i 

3 .11103264-01  .13693754-01 

4 .44373563-01  -.26964517-01 

5 -.73896198-01  -.97876096-01 

6 ,20138093-03  ,802ll493-0i 


PREDICTIVE 

K 

1 


filter  COEFFICIENTS  FOR  PBEST: 

RE  A(KIPBEST)  IM  A(KIPBEST) 
,656078984-00  ,659105584-00 


normalized  CORRELATIONS; 

DELAY  RE  RHO 

0 ,100000004-01 

1 ,656078984-00 

2 -,39805360-02 

3 -.572639414-00 

4 -.747950754-00 

5 -.486176744-00 

6 ,89318492-02 


IM  RHO 
.00000000 
,65910558+00 
,86485063+00 
,56478672+00 
-.68851300-02 
-.49749569+00 
-,64683826+00 
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